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We study the transport of heat in a three dimensional, harmonic crystal of slab 
geometry whose boundaries and the intermediate surfaces are connected to stochastic, 
white noise heat baths at different temperatures. Heat baths at the intermediate surfaces 
are required to fix the initial state of the slab in respect of its surroundings. We allow 
the energy fluxes flown between the intermediate surfaces and the attached baths and 
impose conditions that relate the widths of Gaussian noises of the intermediate baths. 
We show that under those conditions Fourier's law holds in the continuum limit. We 
obtain an exponentially falling temperature profile from high to low temperature end 
of the slab and this nature was already confirmed by Ingen-Hausz's experiment. Profile 
indicates that transport in the steady state involves the processes of conduction and 
radiation of heat. 

Keywords: Heat conduction; Langevin equation 

1. Introduction 

We take a piece of solid bar and establish a steady temperature gradient between its 
two ends. The rest part of the bar remains exposed to the environment. The trans- 
port of heat involves three processes such as, conduction, radiation and absorption 
of heat. In the steady state absorption is absent and the conducted current density 
J(x) obeys Fourier's law: 

J(x) = -kVT(x), (1) 

where VT(x) be the local temperature gradient and k is known as the thermal 
conductivity of the solid. In the steady state the temperature profile of the bar 
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is exponentially falling from high to low temperature end and thi s ve ry nature 
was already confirmed by the experiment performed by Ingen-Hauszl^. Since the 
transport of heat is a non-equilibrium process, it becomes a theoretical challenge 
to derive Fourier's law and the exponentially falling temperature profile from the 
application of non-equilibrium statistical mechanics to a basic model of solid. It is 
observed that the heat transport in one dimension exhibits anomalous^ behaviour. It 
means that the thermal conductivity shows a power law dependence n ~ N a , where 
N be the linear size of the s ystem. Th ere are models which predict divergent (0 < 
a < 1) thermal conductivity 1 1 | 2 | 3 | 4 | 5 | 6 | ^ n one di mens i on Numerical study indicates 
a logarithmic divergence^ of thermal conductivity k ~ In AT in two dimension. A 
power law behaviour^ is also observed in such a system. The validity of Fourier's 
law is confirmed numerically^ in one and two dimensional systems with pining and 
anharmonicity. It has been shown numerically^ that a normal transport takes place 
in the disordered harmonic crystal in three dimension, subjected to an external 
pining potential. The validity of this law is also established using a simulation 
studyl^in three dimensional anharmonic crystal. This study provides a temperature 
profile which is almost linear and very little concave in the upward direction. 

It is worthy to be mentioned of an interesting model of harmonic crystal 
with 'self-consistent' stochastic heat baths^. A rigorous derivation shows that the 
Fourier's law holds in dimension one of this model. However, for a d dimensional 
crystal, when d is very large, the heat conductivity is shown to behave as (Idd) -1 , 
where Id is the coupling to the intermediate reservoirs. Using 'self-consistent' heat 
baths one fixes the flow of energy fluxes between the intermediate sites and the 
attached baths to zero and these conditions give rise to a temperature profile which 
falls linearly from high to low temperature end of the slab. In this paper we use 
the same model in three dimension only, without the on-site binding potential in 
the Hamiltonian. Instead of using 'self-consistent' baths, we allow the heat fluxes 
to be flown between the intermediate surfaces and the attached baths and impose 
conditions among the Gaussian widths of the noises of intermediate heat baths. We 
show that under these conditions Fourier's law holds in the continuum limit and the 
temperature profile exhibits an exponential fall from high to low temperature end of 
the slab. A similar profile was obtained experimentally long before by Ingen Hausz. 
In general, the temperature profile provides an idea about the underlying physical 
processes^ that may take place owing to transport of heat in the steady state limit. 
The transport in our model involves the conduction and radiation of heat. On the 
other hand using 'self-consistent' heat bathi^, one disregards radiation from the 
transport process of heat. 

We organize the paper in the following manner. We introduce the model and 
give the steady state solution of the Langevin's equation in Section 2. We evaluate 
different correlators involving normal coordinates and normal velocities in Section 
3. We give the derivation of the exact temperature profile of the slab and discuss its 
experimental consequences in Section 4. This section also contains the derivation of 
Fourier's law in the continuum limit. The discussion of our results and the conclusion 
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2. Model 

We consider a cubic crystal in three dimension. The form of the Hamiltonian 

ff = £f + £5(*»-*-+«) 2 . (2) 

n n,e 

The displacement field x n is defined on each lattice site n = (n±, 112,^13) where 
ni = 1, • • • ,N, n% = 1, • • • , W 2 , and n 3 = 1, • • • , W 3 . Here e denotes the unit vector 
in the three directions. We choose the value of mass attached to each lattice point 
and the harmonic spring constant as one. We have Langevin's type heat baths that 
are coupled to the surfaces at n\ = 1 and ri\ = N and are maintained at fixed 
temperatures T\ and Tjy (T\ > Tjy) respectively. Apart from these two heat baths, 
the surroundings play a key role to fix the initial state of the slab. Suppose, the 
heat baths are removed from the two end surfaces of the slab. In this situation we 
expect that the temperature of the slab must be equal to the temperature of the 
surroundings in the stationary state limit. To achieve this expectation, in addition 
to two heat baths that are already fixed at the end surfaces of the slab, we require 
to couple N — 2 more Langevin's type heat baths at the intermediate surfaces 
from m = 2 to m = N — 1. To deal with the same problem in one dimension an 
equilibrium phase space distribution has been assumed to fix the initial state of the 
systerrffl Now the equation of motion of a particle at the site n reads 

Xn = -2j(x n - In+e) ~ 7±n + Vn- (3) 
e 

We have chosen the noises to be white and they are uncorrelated at different sites. 
Assuming the process to be Gaussian, the noise correlation is specified by 

M*)TJn'(0> =7z ni 5(t-t')S ajIl ,. (4) 

According to ref. ^^1, z ni is chosen to be proportional to the temperature of the 
corresponding layer. However, we determine them in our model using fluctuation 
dissipation theorem. We use the following periodic boundary conditions for the 
displacement field and the noise: 

£n+(0,lV 2 ,0) W = X n (t) = £n+(0,0,W3)(*) 

r)n+(0,W 2 ,0)(t) = vLt) = ?7n+(0,CW3)W (5) 

These boundary conditions lead to the following expansion for x n (t) and f) a (t): 

fH.it) = -^ m ^f nL ( P ,t)e^, (7) 
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where a be the lattice constant of the crystal, p = (p2,P3) and nj 
substitution of eq. ([6]) and (O into eq. ([3]) we obtain 



JV 



VjkVk - IV] + fj 



(n 2 ,n 3 ). Upon 
(8) 



fc=i 



where the N x N matrix 



V 



(2u% -1 ... ^ 
-1 2cj 2 -1 ... 

-1 2w 2 -1 '•• 



V o 



(9) 



-1 2cjg/ 



and 



^(p) = l + 2sin 2 (^) + 2sin-(i^). 



2/-P3«x 



(10) 



We have assumed here that yo(P7 f ) = = j/at +1 (p, f). To solve eq. (0 we diagonalize 
the matrix V. The solution of the N order polynomial equation V — a 2 l\ — gives 
the eigenvalues of V as 



a£(p) = 2wg(p)+2coB 



7V + 1 



(11) 



The j-th component of the normalized eigenvector corresponding to the eigenvalue 
a\ is given by 



(fc) 



(— 1) J ^ sin ' 



N + 1 



7V + l 



(12) 



The diagonalizing matrix A thus reads as Ajk = a^' such that A T A = I and 
A T VA — a 2 , where (a 2 )jk — a? <5jfc- We introduce a new set of coordinates £j as 



JY 



y.> 



fe=i 



The equation of motion of £j in matrix form can be written as 

£ = -« 2 e - 7^ + /, 



(13) 



(14) 



where / = A T f. Since the matrix V is diagonal with respect to £j (j = 1, • • • , iV), we 
henceforth refer them as normal coordinates. It is evident from this equation that 
each normal mode is acted on by dissipative force and noise offered by the heat 
baths. Moreover, one normal mode is coupled to the other through noise terms. 
The couplings among normal coordinates lead to the flow of energy from one mode 
to the other and thereby leading to a non-ballistic flow of heat in the slab. In the 
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steady state limit (t » I/7) we are interested in the particular solution of the 
equations of motion of £. We use the Fourier transform of 

f^He*"* and fj(t) = / — (15) 

in eq. (fl4|) and obtain 

jP,w = 2 : — ■ ( 16 ) 

lj z — ag(p) — 170; 



Now upon substitution of eq. (fl6| into (fT5|) we obtain 

N pea j i^i 

&*.<> -g/,g M .-^ C p)-^ - 'Afe")- < 17 » 



Now the use of eq. (O and (fTS"]) in ((4]) gives the noise correlation in frequency and 
wave- vector space as 

{fj(p,u)f k (p',u')) =2wiZj6j,k5(L)+u')5 p+p < :Q . (18) 



3. Correlation functions 

For the purpose of calculating physically observable quantities in this model, such 
as temperature profile and heat current density, we require to obtain the correlation 
functions involving normal coordinates and velocities. First we compute the corre- 
lation between normal coordinate and normal velocity using eq. (fl7|) and (fT8|) . We 
perform a frequency integration using delta function in frequency space and obtain 
the correlation as 

N 

& 1 (P,4(p't'))=7^(t-i')W',oE tt f l) «? ) ^ (19) 

3=1 

where 



Ic(t - 



; 2ir (u> 2 - a 2 ki - i^/uj) (uj 2 - a\ 2 + ijuj) ' 
Performing the integration over u> we obtain 



(20) 



IS f ) = J^l!L[l>(t t')9(t f ) + I<(t t')6(t' t)}, (21) 
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where 



A d (/3i,/? 2 ) 


= (cos/?i - cos^ 2 ) 2 + 7 2 (2^ 2 (p) + 


cos /3i 


+ cos /3 2 ), 




= 2(cos/?i — cos fa) cos(ojfe 1 \t — t'\) 




{(4wq + 3cos^i 






+ cos/3 2 )sin(o; fcl |i-t'|)}, 






If{t-t') 


= 2(cos /3i — cos /3 2 ) cos(ojfc 2 \t — t'\) 


2_ 


{(4<Jq +cos/3i 






+3cos/3 2 )sinK. 2 |t-t'|)}, 






/3l,2 


= TTkl, 2 /(N + 1), 








= \R 2 - 7 2 /4- 







(22) 

(23) 

(24) 
(25) 
(26) 

It is clear that I c (t ~ t 1 ) — > 0, when \t — t'\ 00. When t = t' 

We compute the equal time velocity- velocity correlation using eq. (fTTj) and (fT8"|) and 
after performing the algebraic manipulation obtain the correlation as 

,1 . , 7 2 2w§ + cos/3i +cos/3 2 

(^(p,m k2 (-p,t)) = — — — 

N 

x^ iS in(jA)sin(j&). (28) 
4. Temperature profile and Fourier's law 

Consider a particle on the surface at n^. The mean square velocity of the particle 
in the steady state limit reads 



1 N 

= W^rE E a^a^(i kl (p,t)ik 2 (-P,t)), (29) 
p fei,fe 2 =i 

where we have used eq. © and (ITcH) in the last step. Then we substitute eq. (IT21 
and (|28[) in the above equation and evaluate p 2 and ^3 sum in the continuum limit. 
According to this limit, when a — ¥ and W2.3 — > 00, oW2,3 remain fixed and the 
discrete sums over p 2 and p% are converted into integrals: 

isr /?/;*■*•■ « 3o) 

Evaluation of the integrals^l over p 2 and P3 gives 

iV 

^ 2 n (n 1 ) = 2^C„ 1J z, (31) 
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where 

Cm j = nyqrn2 X! a(^i' sin ( ni/3i ) sin ( ni ^ 2 ) sin <jfii) 

^1,^2 — 1 

xsin(j/3 2 ), (32) 
ACSx.jSa) = A(/3!,/3 2 ) - {(cos A - cos/3 2 ) 2 

xF(l/2, 1/2, 1; (4 7 2 /A(/3 1 ,^)) 2 )}, (33) 
= (cos ft - cos/3 2 ) 2 + 7 2 (6 + cos ft + cos ft). (34) 

The N x N matrix C has the following properties: 

CjM = Ckj and Cj t h — CN+i-j,N+i-k- (35) 

To thermalize the layers it is required that the average kinetic energy of a particle 
on a layer in the steady state limit coincides with its energy at thermal equilibrium. 
This requirement gives the following N simultaneous equations of Zk : 

N 1 

X]Cj,fe2fe = -T i , (36) 
fe=i 

where we have chosen the Boltzmann constant — 1. Tj is the temperature of 
the j-th layer. A ccor ding to this equation, if Cj.k — Cj,jSj,k, then only Zj will be 
proportional to tJ^. This equation also indicates that the system in the steady state 
consists of TV different, thermally equilibriated sub-systems which are characterized 
by the temperatures T\, ■ ■ ■ , Tjv. In eq. (|36[) . T± and Tjy are given. Apart from these 
two Tj for 2 < j < N — 1 and 2j for 1 < j < N are remained unknown. We require 
to know additional conditions which will facilitate us to obtain the temperature 
profile of the slab and the ZjS as a solution of these N equations. 

Since T\ > Tjv, there will be a rectilinear flow of heat from the left boundary 
surface at ni = 1 to the right boundary surface at n\ = N of the slab. Heat current 
density j n from the lattice site n to n + §i, where ei = (1, 0, 0), is given bjEl 

in = ^((^n+ej - X n )(i n+&1 + X n )) (37) 

The average heat current density per boncffl] 

N-l W 2 W 3 

W 2 W 3 (N ~ 

n\ — 1 Tl2 — 1713 — 1 

We substitute eq. ([6]) and (1131) in J and after performing the summations over n 2 
and 713 obtain the average heat current density per bond in the steady state limit 
as 

1 N N-l 

J- t V V Vfa (fel) - a {kl) Ma {h2) + a^) 

JV ; p fei,fe>=I m=l 

x^*i(p,*)^(-p,t)>- ( 39 ) 



J = Ld EEEv (38) 
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We now use eq. (IT21 and evaluate the sum over n\. The result is given in the 
following: 
N-x 

E («4Si " < l) )(^ 1 2 ll + < 2) ) = ]vTl(l - (-l) fel+fc2 )sin/3 lS in/3 2 

n i — l 

r 1 

-COS /?2 — COS (3\ 

We use this result and also the eq. (fT9|) and ([27)) in eq. (|39|) . Then we evaluate the 



- 1 



(40) 



discrete sums over p2 and in the continuum limit and obtain 

N 

J = -j^rjJ2 z ^( N ^^ ( 41 ) 

3=1 



where 



1 A (1 - (-l)fci+fe) 



Jj(JV,7) = / jv + 1 N 2 2^ a^i /3 2 ) sm O'/ ? i) sin (j/ 3 2)sin/3isin/3 2 

fci ,k2 — 1 

xirQ,i,l;(4 7 2 /A(^^ 2 )) 2 ). (42) 

Using the property that Jjv+i-j (iV, 7) = —Ij(N,j) and assuming that AT be an 
even number, eq. (|41[) simplifies to 

AT/2 

J = -JVriEfe-^+H)^(^7). (43) 
3=1 

It is knowrff^ that there is also an energy flow between any intermediate surface 
and its attached heat bath. Since we are not using the self consistency condition, 
this intermediate energy fluxes will remain non-zero and affect the per bond average 
current density J. Consequently, J, as is evident from eq. (|4"5)l , is dependent on the 
variables, ^/ 7 Z2, • • • , ^/jzn-i , the Gaussian widths of the noise functions of the 
intermediate heat baths. Assume that the widths are such that 

Z j = z N+ i-j for 2 < j < N/2. (44) 

The above conditions imposed on ZjS are physically interesting firstly because of the 
fact that the use of them leads to a heat current density which obeys Fourier's law in 
the continuum limit. Secondly, the use of them leads to a temperature profile whose 
nature is similar to what has been observed in Ingen-Hausz's experiment. Thirdly, 
these conditions, unlike the self consistency conditions used in Ref. [12l allow the 
radiation to take place in the steady state limit of the transport process. Now the 
use of these conditions in eq. (|36[) gives the following relations: 

1 ( T\ — T/y \ 
Zl ~ ZN = 2 Um-CaJ ' (45) 
Tj T N+1 .j = C ^- CN + l -^ { Tl - T N ) (46) 
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Fig. 1. Temperature profile of the slab for 7 = 0.01. 

for 1 < j < N/2. Moreover, the use of eq. (l44l) in (l36l) reduces N equa- 
tions into following N/2 simultaneous equations of the independent variables Zj 
(j = !,■•• ,N/2): 

N/2 1 / r \ 

E Ci,k*k = 2 \Ti + Cii ^ Cni (Ti ~ T N )J , (47) 
where j = 1, • • • , N/2 and 

Cj.k = Cj : k + C^Ar+i-fc- (48) 

To solve eq. (|4T)1 for Zj (1 < j < N/2), we require to know the temperature profile 
of the slab. The temperatures of the boundary surfaces T± and Tn are known and 
apart from these two, the relations given by eq. (|46[) provide information about the 
temperatures of the intermediate surfaces. The profile which is consistent with those 
relations must be Tj = Tq + Cj t \{T\ — Tjv)/(Ci ) i — Cjv.i), where To is a constant 
independent of j. T will be such that Tj equals to T\ and Tjq for j equals to 1 and 
N respectively. So it gives the temperature profile as 



T 



^ Ci,i + Cn,i Ti — T N 



Tj + T N a 



r r t r > r ^T< — *- Tl ~ Tn ^ (49) 

^l,l-L'JV,l-'l+-'ArJ ^ ^1,1 — Ojv,i 

where 1 < j < N. Moreover, the profile is in conformity with the zeroth law of 
thermodynamics which tells that Tj — T\ for all j when T\ = Tn- The temperature 
profile is plotted in Figure ((TJ). Since the profiles become size independent for N > 
300, the open squares remain coincident with the open circles till j = 300. We use 
eq. (j49"l) in (|47|) to solve Zj for j = ,N/2 and plugging these solutions into 
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eq. (HU) and (14"51 the rest of Zj for j lying in the range N/2 + 1 < j < N are 
obtained. The result of our numerical evaluations are given in the following table. 
The Table[T]indicates that Zj are positive definite and size independent for N > 300. 



Table 1. Zj (1 < j < N) for differ- 
ent N and for 7 = 0.01. 



N 




Zj for 2 < j < TV 


300 


153.61 


1.99 


350 


153.45 


2.0 



We use the conditions of eq. (|44|) into (|43| and then using eq. (|45|) obtain 

J = ~TuTt~ TT ( ^~l N ) Ii(N,l). (50) 



2(N-1) VCi,i-Civ,i, 

I\ (N, 7) is zero if k\ and ki simultaneously take even integer values or odd integer 
values. Assuming that N be an even number and using the fact that the summoned 
of eq. (|42j) is symmetric in respect of the interchange of 0i and /?2, we rewrite the 
double sum of I±(N,j) as 

4 sin 2 B-i sin 2 02 /ll , o , . ,% ~ >*o\ 

w = w+w^-ikjt F (^ l! (47 /A(A ' A)) ) ■ (51) 

where /3i = 2nj 1 /(N + 1) and /3 2 = 7r(2j2 — l)/(-ZV + 1). Again in the continuum 
limit we convert this double sum into integrals. Defining the integration variables 
in this limit as #1,2 = In 31,2 1 {N + 1), we convert the discrete sums into integrals: 

2 N/2 i r 

?£-►-/ d9 ^- ^ 

1 , 7T ,/n 



AT-. . . 

31,2=1 



Ii (N, 7) thus takes the form 
3i(7) = Jim h(N,j) 



1 r r sin 2 61 sin 2 2 /ll , , (A/ , 

<*M <*0 2 a //| n \ ~ F -,-,1; (4 7 2 /A(gi i g 2 )) 2 ■ (53) 



^ 2 ,/ A(0i,0 2 ) V2'2 

Proceeding in the similar manner we obtain the following form of C\.\ — Cn,i in 
the continuum limit as 

32(7) = nm (Ci,i - Cn.i) 

JV->oo 

_ 1 r so r in „:„2 fl -2 a A (°l, 



r2 



f d0! f d02 sin 2 0! sin 2 02^44- (54) 



The steady state current density per bond thus reads in the continuum limit as 

(2\ - 7W) 
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Fig. 2. Plot of k as a function of 7. 



where the conductivity 



k = 



75i(7) 
252(7)' 



(56) 



Here k is found to be finite and independent of the size of the system. So, as a 
consequence of the conditions given by eq. (|44p . Fourier's law holds in the continuum 
limit. The thermal conductivity n, as given by eq. (|56l) . is plotted as a function 
of 7 in Figure [2] Here 7 appears as a constant in the dissipative force term of 
the Langevin's equation. Physically it denotes a viscous force experienced by the 
particles of Brownian like of the slab of crystal owing to collisions with the particles 
of fluid which seems to constitute the heat bath J^ I ^ [ Moreover, owing to collisions, 
particles of the slab will experience random forces or noises. Noise does a significant 
portion of its work in overcoming the viscous drag experienced by a particle of 
the slab and the rest will be flown as heat energy to the neighbouring particles. 
Therefore, the increase of 7 reduces the rate of flow of heat from a particle to its 
neighbours and thereby reducing the thermal conductivity of the system. Hence, it 
justifies reasonably the nature of variation of k with 7 as shown in fig. [2] 

According to eq. (|49| Tj depends linearly on Cj t \. We fit Cj } \ with the following 
exponentially falling function such a way that 



— a o + a i e 



-b\3-l\ 



(57) 



where ao, a\ and b are supposedly N dependent parameters. This parametrized 
form satisfies all the properties of Cj t k given by eq. (|35|) . The numerical values of 
the parameters for different N are given in the following table. The substitution of 
eq. ([57j) into eq. (j49|) gives the temperature profile of the slab in the thermodynamic 
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Table 2. Fitting parameters involved with Cj i for 
different N and for 7 = 0.01. 



N 


a 


01 


b 


500 


1.48 x 10~ b 


4.9868 x lO" 3 


0.03220 


600 


1.48 x 10~ b 


4.9872 x lO" 3 


0.03212 


1000 


7.99 x 10~ b 


4.9881 x W- 3 


0.03193 


2000 


3.71 x 10~ b 


4.9887 x lO" 3 


0.03181 


3000 


2.41 x 10~ b 


4.9889 x 10~ a 


0.03177 



limit as 

Tj =T N + {T 1 -T N )c- b ^ 1 \. (58) 

It is clear that b is the only parameter which appears in the temperature profile of 
the slab in the thermodynamic limit and the Table [2] indicates that it approaches to a 
constant value in this limit. The expression for the profile given by eq. (|58|) is similar 
to what has been obtained in the Ingen-Hausz's experimental. Now it is possible to 
figure out from the temperature profile of the slab the probable physical processes 
that may take place when heat transports through it. According to Refs. [HJ the 
transport of heat in the steady state limit involves the processes of conduction and 
radiation of heat. The situation, when the radiation is absent, the temperature falls 
linearly from left to right boundary surface. On the other hand, if the radiation 
is allowed to take place, temperature falls exponentially and this very nature is 
observed in the Ingen-Hausz's experiment. Thus in this paper, allowing the energy 
fluxes to be flown between the intermediate surfaces and the corresponding heat 
baths and the use of the conditions given by eq. (J44J) , we are basically allowing 
the radiation to take place along with the conduction of heat in the steady state 
limit. If these energy fl uxes are not allowed to be flown then a linear variation of 
temperature is obtained^. It implies that the self consistency conditions do not 
allow the radiation to take place during the steady state transport of heat in the 
slab. 



5. Discussion and conclusion 

According to this model, the additional heat baths coupled at the intermediate 
surfaces and the two baths attached at the boundaries, constitute the surroundings 
of the slab geometry. The additional heat baths are essentially required to fix the 
initial state of the slab. When t » I/7, the slab will attain a unique steady state 
which is a unison of N different, thermally equilibriated sub-systems, characterized 
by the temperatures T\, ■ ■ ■ ,7V- Instead of using self consistent reservoirs^', we 
allow the heat fluxes flown between the intermediate surfaces and the corresponding 
attached heat baths and impose the conditions that the widths of the Gaussian 
noises of j-th and N + 1 — j-th baths are same for j = 1, • • • ,N/2, where N 
being assumed to be even. Those conditions lead to interesting physically admissible 
consequences. We show using them that Fourier's law is satisfied in the continuum 
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limit. Moreover, using those conditions we also obtain a temperature profile, which 
unlike showing a linear variation of Refs. 12 falls exponentially from left to right 
end of the slab and this very nature of the profile is also in conformity with that 
of Ingen-Hausz's experiment. Hence, according to Refs. [TH our model allows the 
conduction and radiation processes to take place during the steady state transport 
of heat in the slab. 
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